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Abstract 

A general chemical kinetics and sensitivity analysis code 
for complex, homogeneous, gas-phase reactions is described. 
The main features of the code, LSENS (the NASA Lewis 
kinetics and sensitivity analysis code), are its flexibility, effi- 
ciency and convenience in treating many different chemical 
reaction models. The models include: static system; steady, 
one-dimensionaj, inviscid flow; incident-shock initiated 
reaction in a shock tube; and a perfectly stirred reactor. In 
addition, equilibrium computations can be performed for 
several assigned states. An implicit numerical integration 
method (LSODE, the Livermore Solver for Ordinary Differ- 
ential Equations), which works efficiently for the extremes 
of very fast and very slow reactions, is used to solve the 
“stiff’ ordinary differential equation systems that arise in 
chemical kinetics. For static reactions, the code uses the 
decoupled direct method to calculate sensitivity coefficients 
of the dependent variables and their temporal derivatives 
with respect to the initial values of dependent variables 
and/or the rate coefficient parameters. Solution methods for 
the equilibrium and post-shock conditions and for perfectly 
stirred reactor problems are either adapted from or based on 
the procedures built into the NASA code CEA (Chemical 
Equilibrium and Applications). 

Introduction 

This paper describes a FORTRAN code, LSENS [1-3], 
the NASA Lewis chemical kinetics and sensitivity analysis 
code, which has been developed for solving complex, homo- 
geneous, gas-phase, chemical kinetics problems. The moti- 
vation for this work is the continuing interest in developing 
detailed chemical reaction mechanisms for complex reac- 
tions such as the combustion of fuels and pollutant formation 
and destruction, and reduced mechanisms for specific appli- 
cations (e.g., computational fluid dynamic simulations). 

Mathematical descriptions of chemical kinetics problems 
constitute sets of coupled, nonlinear, first order ordinary dif- 
ferential equations (ODEs) [4]. The number of ODEs can be 
very large, because of the numerous chemical species 
involved in the reaction mechanism. Further complicating 
the situation are the many simultaneous reactions needed to 
describe the chemical kinetics of practical fuels. For exam- 
ple, the mechanism describing the oxidation of the simplest 
hydrocarbon fuel, methane, involves over 25 species partici- 
pating in nearly 100 elementary reaction steps [5]. 


Validating a chemical reaction mechanism requires repet- 
itive solutions of the governing ODEs for a variety of reac- 
tion conditions. Analytical solutions to the systems of ODEs 
describing chemistry are not possible, except for the simplest 
of cases, which are of little or no practical value. Conse- 
quently, there is a need for fast and reliable numerical solu- 
tion techniques for chemical kinetics problems. 

In addition to solving the ODEs describing chemical 
kinetics, it is often necessary to know what effects variations 
of either initial condition values or chemical reaction mecha- 
nism parameters have on the solution. Such a need arises in 
the development of reaction mechanisms from experimental 
data [6]. The rate coefficients are often not known with great 
precision and in general, the experimental data are not suffi- 
ciently detailed to determine accurately the rate coefficient 
parameters. The development of a reaction mechanism is 
facilitated by a systematic sensitivity analysis, which pro- 
vides the relationships between the predictions of a kinetics 
model and the input parameters of the problem [7-9]. 

The LSENS code is designed for flexibility, convenience 
and computational efficiency. A variety of chemical reaction 
models can be considered. The models include: (1) static 
system; (2) steady, one-dimensional, inviscid flow; (3) reac- 
tion behind an incident shock wave, including boundary 
layer correction; and (4) the perfectly stirred (highly 
back-mixed) reactor. In addition, computations of equilib- 
rium properties can be performed for the following assigned 
states: (1) enthalpy and pressure; (2) temperature and pres- 
sure; (3) internal energy and volume; and (4) temperature 
and volume. For static problems the code computes sensitiv- 
ity coefficients with respect to the initial values of the depen- 
dent variables and/or the three rate coefficient parameters of 
the chemical reactions. A unique feature of the code is the 
capability of solving multiple (essentially limitless) prob- 
lems in a single computer run. 

LSENS is the product of an ongoing effort at NASA 
Glenn to develop and upgrade algorithms and computer 
codes for chemical kinetics applications. The first NASA 
Glenn kinetics code, GCKP (General Chemical Kinetics 
Program [10]), was developed by Bittker and Scullin in 
1972. This code uses the implicit integration method of 
Tyson [11], with some modifications due to Scullin [10]. 
Since then, many integration methods and computer codes 
have been developed for the efficient integration of ODEs in 
general and chemical kinetics equations in particular (see 
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[12] for details). GCKP was replaced with GCKP84 [13], 
which uses the GEAR package [14], as modified by Zeleznik 
[15] to integrate the chemical kinetics ODEs. The work of 
Radhakrishnan [16-19] showed that for chemical kinetics 
applications, LSODE [20,21] was the most efficient and 
accurate code among the several examined. The solver in 
GCKP84 was therefore replaced with LSODE, and the novel 
capability of performing sensitivity analysis for nonisother- 
mal problems (i.e., combustion kinetics) was incorporated 
into the new code, GCKP86 [22], The sensitivity analysis 
computations use the decoupled direct method [1,12,23,24], 
as implemented by Dunker [23] for isothermal kinetics and 
modified by Radhakrishnan [24] for combustion kinetics. 
This method has shown greater efficiency and stability, with 
equal or better accuracy, than other methods of sensitivity 
analysis [1,12,23,24]. Subsequently, many improvements 
and new options were incorporated into GCKP86, and the 
present version, LSENS, replaces the previous NASA gen- 
eral chemical kinetics codes GCKP [10], GCKP84 [13], and 
GCKP86 [22]. 

Description of Code and Capabilities 

The LSENS code has been designed for the following 
reaction models and computations: (1) Static reaction at con- 
stant density or with pressure assigned either as a polynomial 
function of up to third degree in time or in tabular form; (2) 
One-dimensional flow reaction with an assigned pressure or 
area profile, either as (a) a polynomial function of up to third 
degree in time or axial distance or (b) in tabular form as a 
function of time or distance; (3) Static or flow reaction with 
an assigned temperature profile, either as (a) a polynomial 
function of up to third degree in time or distance or (b) in 
tabular form as a function of time or distance; (4) Sensitivity 
analysis for a static reaction; (5) Equilibrium computations 
for the following assigned states: (a) enthalpy and pressure, 
(b) temperature and pressure, (c) internal energy and vol- 
ume, and (d) temperature and volume; (6) Reaction initiated 
by an incident shock wave, including computation of the fro- 
zen and equilibrium post-shock states; and (7) Reaction in a 
perfectly stirred reactor, with either an assigned mass flow 
rate or reactor temperature. This problem may be followed 
by a one-dimensional flow' reaction (i.e., in a plug flow reac- 
tor). 

Any static or one-dimensional flow' kinetic reaction prob- 
lem for which temperature is not assigned may either be adi- 
abatic or have a prescribed rate of heat exchange with its 
environment. For flow problems, the independent variable 
can be either time or axial distance. In addition, the pressure, 
area, and temperature profiles may be specified as functions 
of time or distance, irrespective of the independent variable. 

LSENS was developed on the NASA Lewis Research 
Center's IBM 370/3033 computer using the TSS operating 
system (OS) and the Amdahl 5870 computer using the UTS 
OS. It has also been successfully executed on the following 


computer systems: NASA Lewis Research Center's Amdahl 
5870 using the VM/CMS OS, Cray-X/MP/2/4 using the 
COS and UNICOS operating systems and the CFT and 
CFT77 compilers, Cray-Y/MP/8/6128 using UNICOS 6.0 
and CFT77, the Alliant FX/S, the Convex C220 minicom- 
puter using the Convex 8.0 OS, and the VAX 1 1/750, 1 1/780, 
1 1/785, 6320, 8650, 8800, and 9410; NASA Ames Research 
Center's Cray-2 and Cray-Y/MP, using UNICOS and CFT77 
and Cray-C90, using UNICOS and the F90 compiler; the 
SUN SPARCstation 1 using the Sun 4,1 OS; several IRIS 
workstations using the IRIX 4.0.1 OS and F77 compiler; the 
IBM RISC System/6000 using the AIX 3.1 OS and the XLF 
and F77 compilers; the Pentium I, using MS Windows 95 OS 
and the Digital Visual Fortran 90; and the Pentium II running 
the Linux OS and the Fortran compiler included therein. 

The code consists of a MAIN program, 59 subprograms 
and a BLOCK DATA module. Five of these subprograms are 
function routines; all others are subroutines. The code uses 
the following intrinsic and external routines: ALOG, DABS, 
DBLE, DEXP, DFLOAT, DLOG, DM AX 1 , DMIN1, 
DSIGN, DSQRT, EXP, FLOAT, IABS, IFTX, MAX0, MIN0, 
MOD, READ, SNGL, and WRITE. Namelists and statement 
functions are also used. Finally, the code calls the system 
clock to obtain the total CPU time used since job initiation. 

The different subprograms that comprise the LSENS 
package are arranged in three different groups. The first 
group contains the MAIN program and those related to ther- 
modynamic, transport and kinetics computations. This group 
also includes the numerical procedures for the equilibrium 
and post-incident shock states and perfectly stirred reactor 
(PSR) problems. These solution methods were either 
adapted from or based on the NASA Chemical Equilibrium 
and Applications code CEA [25]. The second group includes 
the subroutines required for sensitivity analysis; several of 
these routines were adapted from the code CHEMDDM 
[26]. The last group contains the subprograms included in 
the code LSODE [20,21], which is used to solve the govern- 
ing ODEs. The BLOCK DATA module is located at the end 
of the code. A detailed flowchart of the MAIN program is 
given in Fig. 1, which essentially illustrates the structure of 
LSENS. The routine KINP (shown in Fig. 1) processes the 
reaction mechanism and much of the input data required to 
solve the problem. The routine SENDDM is the main sub- 
program in that section of the code which performs sensitiv- 
ity analysis. This routine also manages the calls to LSODE. 

LSENS has been arranged as much as possible in a 
“modular ’ fashion, with different subprograms performing 
different tasks. However, to avoid unnecessary work, some 
computations are performed in subprograms other than 
where they naturally belong. Because the code is designed to 
be modular, the number of subprograms is fairly large. How- 
ever, this feature aids in both understanding and, if neces- 
sary, modifying the code. In addition, as improvements are 
made in any calculation procedures or methods built into the 
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code, only the subprograms using these procedures need to 
be replaced. An example is the computation of thermody- 
namic data. The relations built into the code are based on fit- 
ting data over two temperature ranges. Work is now 
underway to extend the temperature range over which the 
calculations are valid [25]. 

Communication among different subprograms is accom- 
plished by means of calling sequences and common blocks, 
which are used extensively in LSENS. The reason for using 
common blocks is to avoid lengthy calling sequences, which 
can significantly deteriorate the efficiency of the program. 
Each subprogram contains type declarations for all variables 
used within it. Such declarations are useful for debugging 
and provide a list of all variables that occur in a routine. The 
type declarations are arranged in a specific order to enhance 
their utility. 

The remainder of this section describes briefly the com- 
putational capabilities of LSENS. Also discussed are conve- 
nience features and calculation procedures used for 
thermodynamic and transport properties. Detailed descrip- 
tions of the code, the theory and numerical solution proce- 
dures built into it, and its usage and accuracy and efficiency 
comparisons with other methods and codes are given by 
Radhakrishnan and Bittker [1-3]. 

Types of .Chemical Reaction and Rate .-Co.effl.cknI 

Many different types of elementary chemical reactions 
are considered. In addition, provision has been made for both 
reversible and irreversible reactions. Each reaction is 
assumed to involve up to a maximum of two different reac- 
tant and two different product species and thus can be writ- 
ten in the general form: 

, , kf ,, 

V { S J + V 2 S 2 Y* V 3 5 3 +v 4 5 4’ 0) 

k r 

where vf is the stoichiometric coefficient (i.e., number of 
moles) of reactant species i in the reaction, V," is the stoichi- 
ometric coefficient of product species i in the reaction, 5/ is 
the chemical symbol for species i, and kj and k r are the for- 
ward and reverse rate coefficients, respectively. In equation 
(1), either species S { or S 4 or both may be either absent or 
the general third-body collisional partner M. Therefore all 
collisional processes, including isomerization and spontane- 
ous activation and deactivation of excited species, are con- 
sidered. In addition, photochemical reactions of the 
following type are allowed: 

hv + v 2 S 2 v 3 S 3 +v 4 S 4 , (2) 


where hv represents a single quantum of radiation absorbed 
by the reactant. 

All reactions are assumed to be elementary, that is, real 
molecular events (e.g., [4]), and so the {v/} and { v/' } are 
integers. Also, all species are assumed to be ideal gases. For 
each reaction, irrespective of its type, the forward rate coeffi- 
cient is usually given by the empirical expression [4] 

k f = Ar ” ex p (-]&)• (3) 

where the preexponential factor A, the temperature exponent 
n and the activation energy E are constants, R is the universal 
gas constant and T is the temperature. For a reversible ele- 
mentary reaction, either the user may specify the reverse rate 
coefficient or it is computed within the code, by using the 
principle of detailed balancing or microscopic reversibility 

[4]: 

k r = k f /K c , (4) 

where K c is the concentration equilibrium constant for the 
reaction [ 1 ]. 

Thermodynamic and Transport Properties 

The thermodynamic properties of the species are com- 
puted by using the empirical equations given by Gordon and 
McBride [25]. For each species i the nondimensional proper- 
ties molar- specific heat at constant pressure ( c p JR ), 
molar- specific enthalpy {h/RT) and the one-atmosphere 
molar- specific entropy (s °j /R) are specified as functions of 
temperature as follows: 


V 

• = a i 1 + a i iT+ a 37 * + a- aT + a ■ ?T , 

(5) 

R 

l 9 i l f /, J I, H It 

h i _ 


( 6 ) 

RT 

’’ 1 2 3 4 5 T 

sT 

11 

. T a i,i T 2 a f.4J . 

, in 7 '+tf, i 2 r+— r + — r +— r +a iT 

(7) 


where a, \ - a x 7 are least-squares coefficients. For each spe- 
cies, two sets of coefficients for use on two adjacent temper- 
ature intervals, 300 to 1000 K and 1000 to 5000 K, are 
included. The data are constrained to give the same results at 
1000 K. The thermodynamic properties of the reacting gas 
mixture are evaluated by applying the Gibbs theorem [27]; 
that is, by simply summing the contributions made by each 
species. 
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The transport properties of the species are computed by 
using the empirical equations given by Svehla f28]. For each 
species i, the dynamic viscosity p, and thermal conductivity 
Ky are specified as functions of temperature, in the form used 
by Maitland and Smith [29]: 

ln M, = + + (8) 

In k,. = cij K ln T+-^+ C -^ + d i K . (9) 

1 T 

The coefficients in these equations are obtained by 
least-squares fitting. The gas mixture dynamic viscosity p 
and thermal conductivity k are computed using formulae 
suggested by Wilke [30,31] and Lindsay and Bromley [32]. 

Heat Transfer Models 

The heat transfer rate between a reacting system and its 
surroundings is, in general, a function of the reacting gas and 
ambient temperatures, as well as flow rate and geometry. It is 
most likely that exact heat exchange rates will not be known 
when modeling an experimental reacting system. Therefore 
the main usefulness of the code will be in determining the 
effects of various assumed heat transfer rates. The default 

method is to prescribe the heat loss rate, Q , (or heat loss rate 
per unit length, Q ) as a polynomial function of the reacting 
mixture temperature: 

Q(otQ') = H Q + H l T+H 2 T 1 + H i T i + H 4 1 A , (10) 

where the { H k } are user-specified constants. For one- 
dimensional flow problems, the code also contains built-in 
procedures for computing Q , by using standard correlations 
for the heat transfer coefficient for laminar and turbulent 
pipe flows (e.g., [33]). 

Chemical Kinetics Problems 

To describe the temporal or spatial evolution of homoge- 
neous chemical reaction systems, equations are needed for 
species concentrations, temperature, density, pressure, and 
possibly, velocity [4]. Mathematical descriptions of combus- 
tion kinetics problems constitute sets of coupled, first-order 
ODEs, which can be generalized as follows: 

■.£<*)' an 

where an underscore represents a vector quantity. In equa- 
tion (11), x is the solution vector with N components, where 
N depends on the problem type, and ^ is the independent 


variable — time, t, or axial distance, x. The solution vector 
contains appropriate mass-specific species mole numbers, 
{<7/ (= number of moles of species i per unit mass of mix- 
ture)}, thermodynamic variables, and when necessary, the 
reacting gas velocity. For clarity in presentation, the depen- 
dence of/on the rate coefficient parameters, the heat transfer 
rate, etc. has been suppressed. 

The initial value problem is to solve for the species 
mass-specific mole numbers, thermodynamic properties and, 
for a flow problem, the velocity at one or more ^ values in a 
prescribed integration interval [^,^e n( iL given the initial con- 
ditions, reaction mechanism and rate coefficient parameters. 

Numerical Integration Procedure 

LSENS uses the double precision version of the pack- 
aged code LSODE [20,21] to solve the ODEs arising in com- 
bustion chemistry. LSODE includes a variable-step, 
variable-order implicit Adams method and a variable-step, 
variable- order backward differentiation formula (BDF) 
method. Both methods are step-by-step methods, and 
advance the numerical solution at each step by 

using linear multistep formulae of the type 

*1 *2 

r„= 'L a jXn- j +h *I.Vjl,-r (12) 

i= i /■ o 

In this equation, the {otj} and {(3y} and the integers 
and K 2 are associated with a particular integration method, 
h n (= - £ w _j) is the stepsize for the current step, Yn-j * s the 

numerical solution vector at and j [= ffXn~j)] 1S the 
approximation to the exact derivative vector (y n .) at 

Extensive experimentation by Radhakrishnan [1] showed 
that the BDF method is consistently superior to the Adams 
method for chemical kinetics applications, and is therefore 
the default method in LSENS. For the BDF method, K x = q 
and K 2 = 0, and equation (12) reduces to 

r„= + ( 13 > 

; = i 

where q is the method order for the step. 

LSODE uses a predictor-corrector scheme, wherein at 

each step an initial guess ) for the solution vec- 

tor at \ n is first produced, and then the guess is improved 

upon by iteration. That is, improved estimates Y ( ” l) 

(m = 1,2,...) are computed until the iteration converges. 
LSODE includes a range of iteration techniques, which can 
be generalized by the recursive relation [21]: 
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(18) 


p [rl m+1) -4 ml ] = X a jl,-j + fc .Po£(C’> - C' ( 14 ) 


(m). 


/("i) 


7= I 


$j,n ~ X a j~ h ”~k + Jl}l ^0 n * 


k = 1 


where the NxY iteration matrix P depends on the iteration 
technique. Radhakrishnan [1] examined the different itera- 
tion techniques built into LSODE for chemical kinetics 
applications, and the most efficient and accurate technique 
was Newton-Raphson iteration, for which 


By using equation (17), equation (18) can be rew ritten as 

PSj. n = £°^j, n -* + Mo^(U- ( 19 > 

k = I j 


P = I-h fJ %J, (15) 

where I is the NxN identity matrix, and the NxN matrix J is 
the Jacobian matrix, with element Jq defined as 

J.j = df/dyj. i,j = 1 N, (16) 

A useful feature of LSODE is that it will estimate the ele- 
ments of the Jacobian matrix by finite-difference approxima- 
tions, if the user chooses not to provide analytical 
expressions for them. However, this method requires (N + 1) 
derivative evaluations for a system of N ODEs and so, as 
shown by Radhakrishnan [1], is much more expensive than 
the use of an analytical Jacobian, especially for large N. 
Therefore, LSENS includes analytical Jacobians for all prob- 
lem types, although the user has the option of using the 
numerical Jacobians generated by LSODE. 


Sensitivity Analysis 


For any static reaction problem the first-order sensitivity 
coefficients {S,y (= 37//3r| ; ‘)} can be computed. Here Y { is the 
numerical solution for the /th (/ = 1 dependent variable 
and r\j is either an initial condition value or a rate coefficient 
parameter (i.e., A , n, E\ see equation (3)). The ODEs for the 
sensitivity coefficients are obtained by differentiating equa- 
tion (11) with respect to T \j and then interchanging the order 
of differentiation with respect to £, and T)j. The result is 


dS . . 3/ 

2sS j = js -j + ^ 


d\ 


(17) 


where J is the Jacobian matrix (eq. (16)) and = 3l73r|;. 


The sensitivity analysis computations use the decoupled 
direct method (DDM) [1,12,23,24]. The DDM solves the 
sensitivity equations separately from, but sequentially with, 
the model equations. The same algorithm that solves the 
ODEs for the chemistry also solves for the sensitivity coeffi- 
cients. The rationale for using the same solver is that the two 
systems of ODEs have the same Jacobian. If the BDF 
method is used to solve equation (17), the resulting formula 
is 


Note that equation (19) can be derived by differentiating 
equation (13) with respect to T y [23]. Although the model 
solution vector X is defined only at discrete points in time or 
space, it can still be considered as a continuous function of 
the {r\j}. Thus the sensitivities {37/3^} calculated from 
equation (19) are the exact sensitivities of X (but not neces- 
sarily of x) with respect to t|y, apart from computer roundoff 
error, provided of course that the Jacobian matrix is accurate. 
Note also that because equation (17) is linear, equation (19) 
is linear. Hence the solution can be obtained explicitly; that 
is, without an iterative predictor-corrector procedure. Thus, 
unlike the calculation procedure for the model solution, there 
is no need to measure or control the error incurred by the 
sensitivity coefficients; of course, one must ensure that the 
model solution is sufficiently accurate, to guarantee accuracy 
of the computed sensitivities. 

Equations (14) and (19) show' the similarity between the 
model and sensitivity equations. The DDM exploits this sim- 
ilarity by alternating the solution of equation (19) with that 
of equation (14). At each step the solution for the 

model problem is advanced. Then the new solution Xi ls 
used in equation (19) to advance the {£j} by the same step. 
The process of advancing X and then the {£j} by the same 
step is repeated until the end of the integration interval. 

An important feature of LSENS is that it can be used to 
generate any number of sensitivity coefficients, from just one 
initial condition or one rate coefficient parameter of one 
reaction to the full set of all TV initial conditions and all 3*NR 
rate coefficient parameters, where NR is the total number of 
reactions. Finally, the linear sensitivity coefficients of the 
temporal derivatives of the dependent variables, that is, 
{3y/3r|^}, may also be computed. Such a capability is use- 
ful, for example, in modeling combustion-acoustics interac- 
tions [34]. 

Provision is made for the user to specify a cutoff level, 
TINY, for the normalized sensitivity coefficients. Any nor- 
malized sensitivity coefficient that is smaller in magnitude 
than TINY is set equal to zero. For rate coefficient parame- 
ters an option to tabulate and print nonzero normalized sensi- 
tivity coefficients in decreasing magnitude is provided. Thus, 
for each dependent variable, the user can obtain a list of 
reaction numbers in order of decreasing importance. 
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Chemic al Equ il ibriu m 

The code has built-in procedures for computing the equi- 
librium composition for the following four assigned states: 
(1) pressure and temperature; (2) pressure and mixture 
mass-specific enthalpy; (3) specific volume and temperature; 
and (4) specific volume and mixture mass-specific internal 
energy. For cases 2 and 4 the equilibrium temperature is also 
determined. The code automatically performs the appropri- 
ate type of equilibrium calculation, depending on the kinet- 
ics problem being solved. 

The calculation procedures for the equilibrium states 
were adapted from the code CEA [25]. The equilibrium state 
is obtained by minimizing either the Gibbs or Helmholtz 
function. In either case, the chemical equilibrium criterion 
that must be satisfied by the system is given by [35] 

NS 

X VjdOj = 0, (20) 

y= i 

where NS is the total number of species in the reacting gas 
mixture and |Xy is here the chemical potential of species j. 
This minimization is subject to several constraints imposed 
by the conservation of atomic species (all cases), the 
assigned enthalpy (case 2), and assigned internal energy 
(case 4). 

The resulting algebraic equations are solved by using a 
descent Newton-Raphson iteration method, which automati- 
cally limits the size of the corrections at each iteration to 
avoid convergence difficulties. Also, to prevent negative con- 
centrations and temperature, the code solves for the loga- 
rithm of the variables. 

Incident Shock 

LSENS includes an option to compute the thermody- 
namic state and velocity behind an incident shock. Two types 
of computations are performed. First, the code solves for the 
“equilibrium” shock conditions, that is, after the shock initi- 
ated reactions have equilibrated. The second calculation pro- 
duces the “frozen” shock conditions immediately after shock 
passage, when the composition is unchanged from its initial 
value. In both cases, the post-shock conditions are obtained 
by solving the mass, momentum and energy conservation 
equations describing steady, inviscid flow of an ideal gas: 


Pi v i = P: v 2' 

(21) 

P l + Plri = P 2 + P2 V 2' 

(22) 

■> 7 

VI v; 


h i + ~2 ~ ft 2 + y- 

(23) 


In these equations, which assume that the coordinate sys- 
tem is attached to the shock, p, K p, and h are respectively 
the mixture density, velocity, pressure and mass-specific 
enthalpy. The subscripts 1 and 2 indicate, respectively, con- 
ditions upstream and downstream of the shock. 

The calculation procedure for both the frozen and equi- 
librium states were adapted from the NASA code CEA [25]. 
Essentially, starting with a guess for the post-shock state, a 
Newton-Raphson iteration procedure, which automatically 
limits the size of the corrections to minimize convergence 
difficulties, is used. To avoid negative variables during the 
solution procedure, the equations are cast in terms of the log- 
arithm of the variables. 

Post-Shock Kingtjg§£lgbkm; Starting with the frozen shock 
state, LSENS follows the progress of the chemical reaction 
in the shocked gas by integrating the ODEs describing 
one- dimensional flow with assigned area, over a prescribed 
time or distance interval. The flow area profile is given by a 
special function, which corrects for frictional losses: correc- 
tions are included for both laminar and turbulent boundary 
layers [2,36-38]. 

Perfectly Stirred Reactor (PSR ) 

Steady state PSR computations can be performed for 
either a specified mass flow rate or a specified reactor tem- 
perature. In the former case, the code solves for the mixture 
composition and temperature at reactor exit. In the latter 
case, the mass flow rate and reactor exit mixture composition 
are computed. The problem type is identified by examining 
the input parameters required for problem solution, and so a 
separate switch need not be set. For steady operation of the 
PSR and both problem types outlined above, the governing 
algebraic equations are obtained from the conservation equa- 
tions for species concentrations and energy: 

/nCOj -o,*) = /= 1, NRS, (24) 

(25) 

In equations (24) and (25), m is the mass flow rate through 
the reactor, W i is the molar rate of production of species i per 
unit volume, V is the reactor volume, NRS is the total num- 
ber of reacting species, £) is the heat loss rate from the reac- 
tor, and the superscript * denotes inlet conditions. The heat 
loss rate can be specified as a polynomial function of up to 
fourth degree in reactor temperature (see equation (10)). 

To solve the system of algebraic equations, LSENS uses 
essentially the same Newton-Raphson iteration technique 
that is built into CEA [25] for computing the equilibrium 
state. In particular, this method automatically limits the size 
of the corrections, to reduce convergence difficulties. Also, 
to avoid negative results, the code solves for the logarithm of 
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the variables. Starting with conditions close to the equilib- 
rium state, a series of perfectly stirred reactor computations 
is performed until the desired mass flow rate or reactor tem- 
perature is reached. This technique is used to minimize the 
possibility of convergence to a false solution — one that is 
mathematically correct but physically unrealistic. The rou- 
tine includes several tests and, when necessary, restarts of 
the calculation to ensure that the solution is physically mean- 
ingful. Finally, it includes tests for possible blowout of the 
chemical reaction within the reactor. 

Combined PSR and plug flow problem: A plug flow calcula- 
tion may be performed after a PSR problem in a single com- 
puter run. Such a two step process is sometimes used as a 
simplified model of a gas turbine combustor. To solve this 
type of problem, the code first performs a PSR calculation 
using the necessary input information for this problem type, 
and then solves the ODEs governing one-dimensional invis- 
cid flow. For the flow problem, the code uses the PSR solu- 
tion as the initial (i.e., inflow) conditions. 

Error Messages 

LSENS contains many error messages — every input 
parameter is tested for legality and consistency with the 
other input variables. For example, the code includes an 
option to check the legality of the reaction mechanism. The 
code verifies that no reaction is duplicated and that each 
reaction satisfies charge and atom balance requirements. 
Reaction duplication may arise because the same reaction is 
written in different forms in different regions of the input 
file. 

If illegal input is discovered, a detailed message is 
printed. Each error message is self-explanatory and com- 
plete. It not only describes the mistake, but also tells the user 
how to fix the problem. During execution, some tests are 
made to ensure that variables are within either given or rea- 
sonable bounds. Any difficulty encountered during execution 
will result in an error exit. A message giving the reason for 
the termination and the name of the subprogram where the 
problem occurred will also be printed. If the computation 
stops prematurely, the user should look for the error message 
near the end of the output file. 

Code Usage 

Two input data files are normally required to execute the 
code. The first one, referred to as the “Standard Thermody- 
namic Data File,” must contain the chemical symbol and 
composition, molar mass and thermodynamic data for each 
species. The second file, the “Problem Data File,” must give 
information about, and data required by, the problem(s) to be 
solved. A third input file, the “Transport Properties Data 
File,” which contains transport property data, is required for 
certain computations. Standard Thermodynamic and Trans- 
port Properties Data Files, containing data from the most 


recent thermodynamic data base of the CEA code [25] for 
many species in the C-H-N-0 system, are supplied with 
LSENS. However, the user has the option of specifying ther- 
modynamic data for some (or all) species, by including them 
in the Problem Data File. Thus the user has a convenient way 
of temporarily changing any thermodynamic data or adding 
new species without modifying the Standard File. However, 
the same thermodynamic data, that is, those given in the 
Standard Thermodynamic Data File and/or the Problem Data 
File, are used for every problem included in the Problem 
Data File. 

The Problem Data File must specify the problem type 
and provide the data needed to solve the problem. This infor- 
mation includes the chemical reaction mechanism, including 
rate coefficient parameters, third-body colli sional efficien- 
cies and inert species, if any. The user must also provide 
details of the problem, including: the independent variable 
for a flow problem; assigned- variable profile; temperature 
data, if a temperature profile is to be specified; heat transfer 
data, if heat transfer is to be calculated; optional information 
such as the input units and output units (cgs, SI, or US cus- 
tomary); and output controls. For a perfectly stirred reactor 
problem, the desired reactor conditions should be specified, 
and data given for the variables that tell the code when inter- 
mediate solutions are needed. The above information must 
be followed by the initial (or inlet) conditions, including 
mixture composition. For static and flow problems, the 
required integration controls must also be specified. Finally, 
if sensitivity analysis is required, the Problem Data File must 
include lists of the parameters with respect to which sensitiv- 
ity coefficients are to be computed and of the dependent vari- 
ables of interest. 

If multiple problems are to be solved during a single exe- 
cution of the code, the required data for the second and sub- 
sequent cases depend on the information, if any, that can be 
used from the previous case(s). The code was designed to 
minimize the amount of input data, especially when they are 
the same for two or more consecutive cases. Thus, for exam- 
ple, the user need not provide the reaction mechanism when 
it either is identical to that used on the previous case or can 
be obtained from it by adding reactions or by modifying rate 
coefficient parameters. Other data for the problem may either 
be included for every new’ case or, in most cases, the data 
most recently specified will be used. 

Illustrative Test Problems 

In addition to the Thermodynamic and Transport Proper- 
ties Files, two Problem Data Files are provided with the 
code. These files serve to illustrate the problem types that 
can be solved by LSENS and the options built into it. The 
data files also help elucidate the construction of the problem 
data file required to execute LSENS. To demonstrate the 
capability of solving multiple problems in a single run, both 
problem files contain several test cases, as discussed below. 
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The first data file contains 16 kinetics-only (i.e., no sensi- 
tivity analysis) test cases, which are summarized in Table 1 . 
In this table NR is the total number of elementary reactions 
in the reaction mechanism, NS is the total number of species, 
and N is the total number of ODEs solved. It is clear from 
this table that the test cases cover a variety of reaction mod- 
els and problem types. The problem size varies from a sim- 
ple case involving one reaction among three species (case 1) 
to one involving 42 species participating in 143 reactions 
(case 1 1). On a Silicon Graphics Indigo 2 computer with the 
R4000 processor chip, LSENS required approximately 14 s 
CPU time to solve all 16 problems. The execution time alone 
(i.e., not including the CPU time for preprocessing and 
input/output) was less than 10 s. On a Pentium II the total 
CPU time (i.e., including preprocessing and input/output) for 
the 16 cases was comparable (approximately 9 s). Additional 
run time information and detailed comparisons with other 
methods and codes is available in Radhakrishnan [1]. 

The second data file contains nine sensitivity analysis test 
cases, which are summarized in Table 2. Here, N P is the total 
number of sensitivity parameters (i.e., initial conditions and 
rate coefficient parameters). The first six cases are isothermal 
problems, for which sensitivity analysis results have been 
published in the literature. Results for cases 1 and 2 were 
also computed analytically. Excellent agreement was 
obtained for all six cases between the LSENS and literature 
results [1]. The last three cases are nonisothermal problems, 
which illustrate the application of sensitivity analysis to 
combustion kinetics. Case 7 describes a simple nonisother- 
mal problem for which the analytical solution is known 
[1,24]. The results produced by LSENS showed excellent 
agreement with the analytical solution for the model problem 
and for the sensitivity coefficients with respect to all four ini- 
tial conditions and the three rate coefficient parameters of the 
reaction. For cases 8 and 9, sensitivity coefficients were gen- 
erated by using finite difference approximations; agreement 
with LSENS was again excellent [1]. 

On a Silicon Graphics Indigo 2 computer with the R4000 
processor chip, LSENS required approximately 10 s CPU 
time for all nine problems. The execution time alone (i.e., 
not including the CPU time for preprocessing and input/out- 
put) was approximately 8 s. On the Pentium II the total CPU 
time (i.e., including preprocessing and input/output) for the 
nine cases was approximately 7 s. As for the kinetics-only 
test cases, the sensitivity analysis method built into LSENS 
has been compared with other methods and computer codes 
[1,12]. For example, case 9, which describes the ignition and 
subsequent combustion of a benzene-oxygen-argon mixture, 
consists of 120 reversible reactions among 39 reacting spe- 
cies and the inert species argon. To solve for the dependent 
variables required 6.3 s of execution time on the Amdahl 
computer, using the UTS operating system and the Fujitsu 77 
compiler. To solve for the dependent variables and sensitivity 


coefficients with respect to all 42 initial conditions (40 spe- 
cies mass-specific mole numbers, density, and temperature) 
and all 360 rate coefficient parameters required 246 s. This 
execution time may be compared to the more than 2500 s 
that would be required by a finite difference method wherein 
parameters are varied one at a time. 

Concluding Remarks 

The NASA Lewis kinetics and sensitivity analysis code, 
LSENS, w'hich has been developed for complex, homoge- 
neous, gas-phase reactions, was described. The code has 
been designed for a variety of reaction models: static system; 
steady, one-dimensional, inviscid flow; reaction initiated by 
an incident shock; and a perfectly stirred reactor. The differ- 
ent capabilities and computational procedures built into the 
code and its usage were briefly described. Standard thermo- 
dynamic and transport properties data are provided with the 
code, as are two problem data files that illustrate the different 
options built into LSENS. 

Details regarding code availability and procurement can 
be obtained from the National Technology Transfer Center, 
Wheeling Jesuit University, 316 Washington Avenue, 
Wheeling, WV 26003. Telephone: (800) 678-6882. (NTTC 
supersedes COSMIC, University of Georgia, in this activity.) 
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TABLE 1 . Description of kinetics-only test problems supplied with LSENS 


Test 

problem 

Description of problem 

Problem 

parameters 



NR 

NS 

N 

1 

Dissociation of bromine in a shock tube. 

I 

3 

5 

2 

Hydrogen-air combustion in supersonic flow with heat transfer. 

36 

18 

19 

3 

Same problem as case 2, but with a larger reaction mechanism. 

37 

18 

19 

4 

Hydrogen-oxygen combustion in subsonic flow with heat transfer. 

18 

8 

11 

5 

Static, rich methane-air combustion at assigned pressure with heat transfer. 

133 

39 

41 

6 

Lean methane-air adiabatic combustion in supersonic flow at constant pressure. 

133 

39 

41 

7 

Same problem as case 6, except that the area obtained from case 6 is assigned as input. 

133 

39 

42 

8 

Same problem as case 6, except that the area and temperature obtained from case 6 are 
assigned as input. 

133 

39 

41 

9 

Constant-density, static methanol-air adiabatic combustion. 

133 

39 

40 

10 

Constant-temperature, constant-pressure, methane-air reaction in supersonic flow. 

143 

42 

42 

11 

Perfectly stirred reactor combustion of a rich propane-air mixture, followed by supersonic 
expansion of the combustion products through a diverging nozzle with heat loss. 

136 

42 

44 

12 

High-temperature, adiabatic air ionization reaction in constant-area subsonic flow 7 . 

12 

12 

15 

13 

High-temperature, high-pressure, static reaction of carbon monoxide-hydrogen mixture at 
constant density and constant temperature. 

23 

13 

13 

14 

Constant-density, adiabatic, static problem involving photolytic ignition of hydrogen-oxygen 
mixture at low initial temperature. 

20 

8 

9 

15 

Methane-air combustion in supersonic flow 7 , with assigned duct area and output required at 
axial locations corresponding to specified values of duct area. 

133 

39 

42 

16 

Same problem as case 15, except here output is required at specified axial positions. 

133 

39 

42 


TABLE 2. Description of sensitivity test problems supplied with LSENS 


Test 

Description of problem 

Problem parameters 

problem 


NR 

NS 

N 

N P 

1 

Constant-volume, isothermal reaction, A = B. 

2 

2 

2 

4 

2 

Constant- volume, isothermal reaction, A = B = C. 

4 

3 

3 

4 

3 

Constant- volume, isothermal pyrolysis of ethane, using a simplified mechanism. 

5 

m 

8 

4 

4 

Constant-volume, isothermal reaction of a methane-oxygen-argon mixture con- 
taining trace concentrations of carbon dioxide and hydrogen. 

44 

■ 

13 

44 

5 

Constant- volume, isothermal oxidation of a formaldehyde -carbon monoxide mix- 
ture. 

25 

15 

14 

25 

6 

Constant- volume, isothermal reaction of a wet carbon monoxide-oxygen-nitrogen 
mixture. 

52 

12 

11 

58 

7 

Adiabatic, constant-pressure isomerization reaction with simplified rate coeffi- 
cient expression that permits analytical solution. 

1 

2 

4 

7 

8 

Constant-pressure, adiabatic, static ignition of a stoichiometric hydrogen-air mix- 
ture, seeded with 0.45% nitric oxide. 

40 

19 

20 

20 

9 

Adiabatic, constant-density, static combustion of a shock-heated, near- 
stoichiometric benzene-oxygen-argon mixture. 

120 

40 

42 

51 
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